#include "spline.h"
#include <QMap>
Spline::Spline()
{

}

Spline::~Spline()
{
}

void Spline::Init(const QVector<double> & vx_in,const QVector<double> & vy_in)
{
    m_nSize = (int)vx_in.size();
    if(m_nSize==0) return;
    m_vx = vx_in;
    m_vy = vy_in;
    SortData();
    SetBoundCondtion();
    CalDiff();
    CalDiaFactor();

}

double Spline::GetInterPol(double x_in)
{
    if(m_nSize==0) return 0;
    int k=0;
    int N = m_nSize -1;
    while(x_in>m_vx[k++]&&k<=N);
    k=k-2;
    if(k<0)k=0;
    double h = m_vDx[k];
    double xa = x_in - m_vx[k];
    double xb = x_in - m_vx[k+1];
    double P1,P2,P3,P4;
    P1=((h+2 * xa) * xb * xb * m_vy[k]) / pow(h,3);			// Hermite插值多项式基函数α0
    P2=((h-2 * xb) * xa * xa * m_vy[k+1]) / pow(h,3);		// Hermite插值多项式基函数α1
    P3=(xa * xb * xb * m_vDf[k]) / pow(h,2);					// Hermite插值多项式基函数β0
    P4=(xb * xa * xa * m_vDf[k+1]) / pow(h,2);				// Hermite插值多项式基函数β1
    double y = P1 + P2 + P3 + P4;
    return y;
}
void Spline::SortData()
{

    QVector<double>  vx;
    //先将vx进行过滤，再进行从小到大排序
    QMap<double, double> vxmap;
    //先将重复的记录过滤掉
    for (int i = 0;i<m_nSize;i++)
    {
        vxmap.insert(m_vx[i],m_vy[i]);
    }
    //再将m_vx的过滤值放进vector里进行排序
    QMap<double, double>::iterator iter = vxmap.begin();
    while(iter!=vxmap.end())
    {
        vx.push_back(iter.key());
        ++iter;
    }
    Fastsort(vx, 0, (int)vx.size()-1);
    m_vx = vx;
    QVector<double>::iterator vxIter = m_vx.begin();
    m_vy.clear();
    while(vxIter!=m_vx.end())
    {
        m_vy.push_back(vxmap[*vxIter]);
        ++vxIter;
    }
    m_nSize = (int)m_vx.size(); //重新调整m_nSize的大小
}
void Spline::Fastsort(QVector<double> &array, int begin, int end)
{
    if (begin < 0 || end < 0 || begin > end)
        return;
    int left = begin, right = end;
    double temp = array[left];
    while (right != left)
    {
        while (temp < array[right] && right > left)
            right--;
        if (right > left)
        {
            array[left] = array[right];
            left++;
        }
        while (temp > array[left] && right > left)
            left++;
        if (right > left)
        {
            array[right] = array[left];
            right--;
        }
    }
    array[right] = temp;
    Fastsort(array, right + 1, end);
    Fastsort(array, begin, right - 1);
}


void Spline::CalDiff()
{
    m_vDx.resize(m_nSize-1);
    for (int i =0;i<m_nSize-1;i++)
        m_vDx[i] = m_vx[i+1] - m_vx[i];
}

void Spline::SetBoundCondtion()
{
    m_dDy0 = 0;
    m_dDyn = 0;
}

void Spline::CalDiaFactor()
{
    QVector<double>      vA(m_nSize);		// 对角线1
    QVector<double>		vC(m_nSize);			// 对角线3
    QVector<double>		vRval(m_nSize);    // 右侧值
    int N = m_nSize -1;
    vC[0]=1;
    vA[N]=1;
    vRval[0] = 3 * (m_vy[1]-m_vy[0]) / m_vDx[0]     - m_dDy0 * m_vDx[0]/2;
    vRval[N] = 3 * (m_vy[N]-m_vy[N-1]) / m_vDx[N-1] + m_dDyn*m_vDx[N-1]/2;

    for(int k=1;k<N;k++)
    {
        vA[k]=m_vDx[k]/(m_vDx[k]+m_vDx[k-1]);	//第一条对角线
        vC[k]=1-vA[k];											//第三条对角线
        //右端项
        vRval[k]=3*(vC[k] * (m_vy[k+1]-m_vy[k]) / m_vDx[k]
        +vA[k] * (m_vy[k]-m_vy[k-1]) / m_vDx[k-1]);
    }
    Cal1stDun(vA,vC,vRval);

}

void Spline::Cal1stDun(QVector<double> &vA_in,QVector<double> &vC_in,QVector<double> &vR_in)
{
    QVector<double>		vB(m_nSize);								// 对角线2
    m_vDf.resize(m_nSize);
    vC_in[0]/=2;
    vR_in[0]/=2;
    int N = m_nSize -1;
    for(int i=1;i<N;i++)
    {
        vB[i]= 2 - vA_in[i] * vC_in[i-1];									// 求第二条对角线
        vC_in[i]/=vB[i];															// 求第三条对角线
        vR_in[i]=(vR_in[i]- vA_in[i] * vR_in[i-1]) / vB[i];           //求y

    }

    vB[N]=2-vA_in[N]*vC_in[N-1];
    m_vDf[N]=(vR_in[N]-vA_in[N]*vR_in[N-1])/vB[N];
    for(int i=N-1;i>=0;i--)
        m_vDf[i]=vR_in[i]-vC_in[i]*m_vDf[i+1];						//求x
}
